% This function computes the asymptotic variance for Theta1 and the
% factors in the SR-approach
function out = getAVarTheta1andFactors(theta1,setup,epsValue)

% Evaluating the filter at optimium
theta1Values        = struc2values(theta1,fieldnames(theta1));
[Q,output,errorMes] = objectFunc(theta1Values,setup);

% Ensuring no transformation is used on the factors in modelFit
output.model.factorSignResOn = 0;

% Setting the dimensions for the pooled data
nx               = size(output.factors,1);
numTheta1        = size(setup.selectTheta1,1);
T                = size(output.factors,2);
nyIndex          = zeros(T,1);
for t=1:T
    nyIndex(t,1) = sum(~isnan(setup.yields(t,:)));
end
N                = sum(nyIndex,1);
residualsPooled  = zeros(N,1);
RvHat            = zeros(T,1);
RvHat_se         = zeros(T,1);
dg_dTheta1       = zeros(max(nyIndex),numTheta1,T);
dg_dx            = zeros(max(nyIndex),nx,T);
dx_dTheta1       = zeros(nx,numTheta1,T);

% Computing dx_dTheta1_tmp
epsP = ones(numTheta1,1)*epsValue;
epsM = ones(numTheta1,1)*epsValue;
if setup.AVarTheta1On == 1
    for i=1:numTheta1
        % A positive change in theta1Values
        errorMes = 1;
        while errorMes == 1
            theta1_p_eps = theta1Values;
            theta1_p_eps(i,1) = theta1_p_eps(i,1) + epsP(i,1);
            [Q,tmp,errorMes] = objectFunc(theta1_p_eps,setup);
            if errorMes == 1
                disp(['Positive step: downscaling epsilon at position ',num2str(i)])
                epsP(i,1) = epsP(i,1)/10;
            end
        end
        output_p_eps(i) = tmp;
        output_p_eps(i).model.factorSignResOn = output.model.factorSignResOn;
        
        % A negative change in theta1Values
        errorMes = 1;
        while errorMes == 1
            theta1_m_eps = theta1Values;
            theta1_m_eps(i,1) = theta1_m_eps(i,1) - epsM(i,1);
            [Q,tmp,errorMes] = objectFunc(theta1_m_eps,setup);
            if errorMes == 1
                disp(['Negative step: downscaling epsilon at position ',num2str(i)])
                epsM(i,1) = epsM(i,1)/10;
            end
        end
        output_m_eps(i) = tmp;
        output_m_eps(i).model.factorSignResOn = output.model.factorSignResOn;
        
        % The gradient
        dx_dTheta1(1:nx,i,1:T) = (output_p_eps(i).factors - output_m_eps(i).factors)/(epsM(i,1)+epsP(i,1));
    end
end
% Auxiliary terms to compute B
nyCount = 0;
autoCov1Residuals = zeros(max(nyIndex),max(nyIndex));
varResiduals      = zeros(max(nyIndex),max(nyIndex));
countFull         = 0;
for t=1:T
    % The value of dg/dx - for fixed Theta1
    select = ~isnan(output.residuals(:,t));
    epsValueX = 1D-8*0+epsValue;
    tmp_p_eps = zeros(max(nyIndex),nx);
    tmp_m_eps = zeros(max(nyIndex),nx);
    for i=1:nx
        factorsGrad      = output.factors(:,t);
        factorsGrad(i,1) = factorsGrad(i,1) + epsValueX;
        scoreFullTmp     = modelFit(output.model,factorsGrad,t);
        tmp_p_eps(:,i)   = scoreFullTmp(1:max(nyIndex));
        
        factorsGrad      = output.factors(:,t);
        factorsGrad(i,1) = factorsGrad(i,1) - epsValueX;
        scoreFullTmp     = modelFit(output.model,factorsGrad,t);
        tmp_m_eps(:,i)   = scoreFullTmp(1:max(nyIndex));
    end
    dg_dx(select,1:nx,t) = (tmp_p_eps(select,:) - tmp_m_eps(select,:))/(2*epsValueX);

    % The value of dg/dTheta1  - for fixed factor values
    if setup.AVarTheta1On == 1
        for i=1:numTheta1
            tmp = (modelFit(output_p_eps(i).model,output.factors(:,t),t)...
                -modelFit(output_m_eps(i).model,output.factors(:,t),t))/(epsM(i,1)+epsP(i,1));
            dg_dTheta1(select,i,t) = tmp(select);
        end
    end
    
    % The estimated values of the noise yield in the homoskedastic case
    RvHat(t,1)    = sum(output.residuals(select,t).^2,1)/(nyIndex(t,1)-nx);
    RvHat_se(t,1) = 1/nyIndex(t,1)*(1/(nyIndex(t,1)-nx)*sum(output.residuals(select,t).^4,1)-RvHat(t,1)^2);

    % The pooled residuals
    residualsPooled(nyCount+1:nyCount+nyIndex(t,1),1)  = output.residuals(select,t);
    
    % co-variance matrix and 1 lag auto-covariance in residuals
    % Estimated starts when all observations are available.
    if t > 1
        if nyIndex(t,1) == max(nyIndex) && nyIndex(t-1,1) == max(nyIndex)
            varResiduals = varResiduals  + output.residuals(:,t)*output.residuals(:,t)';
            
            autoCov1Residuals = autoCov1Residuals + output.residuals(:,t)*output.residuals(:,t-1)';
            countFull = countFull + 1;
        end
    end
    
    % Update of nyCount and selectLag1
    nyCount = nyCount + nyIndex(t,1);
end
autoCov1Residuals = autoCov1Residuals/countFull;
varResiduals      = varResiduals/countFull;

% The estimate of Var(Rv)
RvHatPooled     = sum(residualsPooled(:,1).^2)/(N-T*nx-numTheta1);
RvHatPooled_Var = 1/N*(1/(N-T*nx-numTheta1)*sum(residualsPooled(:,1).^4,1)-RvHatPooled^2);
RvHatPooled_se  = sqrt(RvHatPooled_Var);

% The variance of the score function B and the Hessian matrix A
% Note both versions of B accounts for cross-sectional correlation and
% autocorrelation. By setting wD = 0 and wT = 0, we obtain estimates that
% do not account for any correlation.
B_Theta1Robust     = zeros(numTheta1,numTheta1);    
B_Theta1Homo       = zeros(numTheta1,numTheta1);    
B_Theta1HomoPooled = zeros(numTheta1,numTheta1);     
A_Theta1           = zeros(numTheta1,numTheta1);
PSI                = zeros(numTheta1,T,max(nyIndex));
if setup.AVarTheta1On == 1
    for t=1:T
        for j=1:nyIndex(t,1)
            PSI(1:numTheta1,t,j) = dx_dTheta1(1:nx,1:numTheta1,t)' * reshape(dg_dx(j,1:nx,t),nx,1) + ...
                reshape(dg_dTheta1(j,1:numTheta1,t),numTheta1,1);
        end
    end
end
idx = 0;
if setup.AVarTheta1On == 1
    % We use the symmetry in the cross section dimension => more observations are included in the estimation
    for t=1+setup.wT:T-setup.wT
        for j=1:nyIndex(t,1)-setup.wD
            idx = idx + 1;
            A_Theta1 = A_Theta1 + PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t,j)';
            for kT=-setup.wT:setup.wT
                for kD=0:setup.wD
                    if kD == 0
                        tmpRobust = (1-abs(kT)/(1+setup.wT))*(1-abs(kD)/(1+setup.wD))*...
                            output.residuals(j,t)*output.residuals(j+kD,t+kT)*PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t+kT,j+kD)';
                    else
                        tmpRobust = (1-abs(kT)/(1+setup.wT))*(1-abs(kD)/(1+setup.wD))*...
                            ( output.residuals(j,t)*output.residuals(j+kD,t+kT)*PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t+kT,j+kD)'...
                            +output.residuals(j,t)*output.residuals(j+kD,t+kT)*PSI(1:numTheta1,t+kT,j+kD)*PSI(1:numTheta1,t,j)');
                    end
                    if ~isnan(tmpRobust)
                        B_Theta1Robust = B_Theta1Robust + tmpRobust;
                    end
                    if kT == 0 && kD == 0
                        tmpHomo       = RvHat(t,1)*PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t,j)';
                        tmpHomoPooled = RvHatPooled*PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t,j)';
                    else
                        tmpHomo = (1-abs(kT)/(1+setup.wT))*(1-abs(kD)/(1+setup.wD))*...
                            (output.residuals(j,t)*output.residuals(j+kD,t+kT)*PSI(1:numTheta1,t,j)*PSI(1:numTheta1,t+kT,j+kD)'...
                            +output.residuals(j,t)*output.residuals(j+kD,t+kT)*PSI(1:numTheta1,t+kT,j+kD)*PSI(1:numTheta1,t,j)');
                        tmpHomoPooled = tmpHomo;
                    end
                    if ~isnan(tmpHomo)
                        B_Theta1Homo       = B_Theta1Homo + tmpHomo;
                        B_Theta1HomoPooled = B_Theta1HomoPooled + tmpHomoPooled;
                    end
                end
            end
        end
    end
    B_Theta1Robust     = B_Theta1Robust/idx;
    B_Theta1Homo       = B_Theta1Homo/idx;
    B_Theta1HomoPooled = B_Theta1HomoPooled/idx;
    A_Theta1           = A_Theta1/idx;
    invA_Theta1        = A_Theta1\eye(numTheta1);
else
    B_Theta1Robust = B_Theta1Robust*NaN;
    B_Theta1Homo   = B_Theta1Homo*NaN;
    A_Theta1       = A_Theta1*NaN;
    invA_Theta1    = A_Theta1*NaN;
end
VarRobustTheta1      = invA_Theta1*B_Theta1Robust*invA_Theta1/N;
VarHomoTheta1        = invA_Theta1*B_Theta1Homo*invA_Theta1/N;
VarHomoPooledTheta1  = invA_Theta1*B_Theta1HomoPooled*invA_Theta1/N;   % if sigma_t = sigma for all t
% Note for the iid case, then  VarHomoPooledTheta1 = invA_Theta1*RvHatPooled/N

% Calculating VarFactors
VarFactorsRobust     = zeros(nx,nx,T);
VarFactorsHomo       = zeros(nx,nx,T);
VarFactorsHomoPooled = zeros(nx,nx,T);
autoCov10Factors     = zeros(nx,nx,T);
autoCov01Factors     = zeros(nx,nx,T);
if setup.AVarAdjThetaOn == 1 && setup.AVarTheta1On == 1
    % Adjustment for uncertainty in Theta
    for t=1:T
        % The value of A
        A = dg_dx(1:nyIndex(t,1),1:nx,t)'*dg_dx(1:nyIndex(t,1),1:nx,t)/nyIndex(t,1);
        
        % The robust case
        idx    = 0;
        B_Robust = zeros(nx,nx);
        for j=1:nyIndex(t,1)-setup.wD
            idx = idx + 1;
            % No symmetry in the cross-section dimension
            %for kD = -setup.wD:setup.wD
            %    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t)*output.residuals(j,t)*output.residuals(j+kD,t);
            %end
            % Symmetry in the cross-section dimension
            for kD = 0:setup.wD
                if kD == 0
                    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)*dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t);
                else
                    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)...
                        *(dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t) + dg_dx(j+kD,1:nx,t)'*dg_dx(j,1:nx,t));
                end
            end
        end
        B_Robust = B_Robust / idx;
        
        % The value of B (home case)
        B_Homo = RvHat(t,1)*A;
        B_HomoPooled = RvHatPooled*A;
        
        % The value of D
        D = zeros(nx,numTheta1);
        for j=1:nyIndex(t,1)
            D = D + reshape(dg_dx(j,1:nx,t),nx,1)*reshape(PSI(1:numTheta1,t,j),1,numTheta1);
        end
        D = D/nyIndex(t,1);
        
        % The value of X_theta1_factors (Robust and Homo case)
        C_theta1_factors_Robust = zeros(numTheta1,nx);
        for j=setup.wD+1:nyIndex(t,1)-setup.wD
            for kD=0:setup.wD
                % Use of the symmetry
                if kD == 0
                    C_theta1_factors_Robust = C_theta1_factors_Robust + ...
                        (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)...
                        *reshape(PSI(1:numTheta1,t,j),numTheta1,1)*reshape(dg_dx(j+kD,1:nx,t),1,nx);
                else
                    C_theta1_factors_Robust = C_theta1_factors_Robust + ...
                        (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)...
                        *(reshape(PSI(1:numTheta1,t,j),numTheta1,1)*reshape(dg_dx(j+kD,1:nx,t),1,nx) + ...
                        reshape(PSI(1:numTheta1,t,j+kD),numTheta1,1)*reshape(dg_dx(j,1:nx,t),1,nx));
                end
            end
        end
        C_theta1_factors_Robust = C_theta1_factors_Robust/(nyIndex(t,1)-setup.wD);
        C_theta1_factors_Homo   = RvHat(t,1)/nyIndex(t,1)*reshape(PSI(1:numTheta1,t,1:nyIndex(t,1)),numTheta1,nyIndex(t,1))...
            *reshape(dg_dx(1:nyIndex(t,1),1:nx,t),nyIndex(t,1),nx);
        C_theta1_factors_HomoPooled  = RvHatPooled/nyIndex(t,1)*reshape(PSI(1:numTheta1,t,1:nyIndex(t,1)),numTheta1,nyIndex(t,1))...
            *reshape(dg_dx(1:nyIndex(t,1),1:nx,t),nyIndex(t,1),nx);
        
        % The distribution for the estimation in factors
        c = nyIndex(t,1)/N;
        invA = A\eye(nx);
        VarFactorsRobust(:,:,t) = 1/nyIndex(t,1)*invA*(B_Robust + c*D*VarRobustTheta1*D' ...
            -sqrt(c)*(C_theta1_factors_Robust'*invA_Theta1*D' + D*invA_Theta1*C_theta1_factors_Robust))*invA;
        VarFactorsHomo(:,:,t)   = 1/nyIndex(t,1)*invA*(B_Homo + c*D*VarHomoTheta1*D' ...
            -sqrt(c)*(C_theta1_factors_Homo'*invA_Theta1*D' + D*invA_Theta1*C_theta1_factors_Homo))*invA;
        VarFactorsHomoPooled(:,:,t) = 1/nyIndex(t,1)*invA*(B_HomoPooled + c*D*VarHomoPooledTheta1*D' ...
            -sqrt(c)*(C_theta1_factors_HomoPooled'*invA_Theta1*D' + D*invA_Theta1*C_theta1_factors_HomoPooled))*invA;
        
        % The autocovariance matrix
        if t > 1
            select     = ~isnan(output.residuals(:,t));
            selectLag1 = ~isnan(output.residuals(:,t-1));
            % cov(u_t+1,u_t)
            autoCov10Factors(:,:,t) = 1/nyIndex(t,1)*1/nyIndex(t-1,1)*invA*dg_dx(1:nyIndex(t,1),1:nx,t)'...
                *autoCov1Residuals(select,selectLag1)*dg_dx(1:nyIndex(t-1,1),1:nx,t-1)*invALag1;
            % cov(u_t,u_t+1)
            autoCov01Factors(:,:,t) = 1/nyIndex(t,1)*1/nyIndex(t-1,1)*invALag1*dg_dx(1:nyIndex(t-1,1),1:nx,t-1)'...
                *autoCov1Residuals(selectLag1,select)*dg_dx(1:nyIndex(t,1),1:nx,t)*invA;
        end
        
        % Update of nyCount and invA
        invALag1 = invA;
    end
else
    % No adjustment for uncertainty in Theta
    for t=1:T
        % The value of A
        A = dg_dx(1:nyIndex(t,1),1:nx,t)'*dg_dx(1:nyIndex(t,1),1:nx,t)/nyIndex(t,1);
        
        % The robust case
        idx    = 0;
        B_Robust = zeros(nx,nx);
        for j=1:nyIndex(t,1)-setup.wD
            idx = idx + 1;
            % No symmetry in the cross-section dimension
            %for kD = -setup.wD:setup.wD
            %    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t)*output.residuals(j,t)*output.residuals(j+kD,t);
            %end
            % Symmetry in the cross-section dimension
            for kD = 0:setup.wD
                if kD == 0
                    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)*dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t);
                else
                    B_Robust = B_Robust + (1 - abs(kD)/(1+setup.wD))*output.residuals(j,t)*output.residuals(j+kD,t)...
                        *(dg_dx(j,1:nx,t)'*dg_dx(j+kD,1:nx,t) + dg_dx(j+kD,1:nx,t)'*dg_dx(j,1:nx,t));
                end
            end
        end
        B_Robust = B_Robust / idx;
        
        % The value of B (home case)
        B_Homo = RvHat(t,1)*A;
        B_HomoPooled = RvHatPooled*A;
        
        % The distribution for the estimation in factors
        invA = A\eye(nx);
        VarFactorsRobust(:,:,t)     = 1/nyIndex(t,1)*invA*B_Robust*invA;
        VarFactorsHomo(:,:,t)       = 1/nyIndex(t,1)*invA*B_Homo*invA;
        VarFactorsHomoPooled(:,:,t) = 1/nyIndex(t,1)*invA*B_HomoPooled*invA;
        
        % The autocovariance matrix
        if t > 1
            select     = ~isnan(output.residuals(:,t));
            selectLag1 = ~isnan(output.residuals(:,t-1));
            % cov(u_t+1,u_t)
            autoCov10Factors(:,:,t) = 1/nyIndex(t,1)*1/nyIndex(t-1,1)*invA*dg_dx(1:nyIndex(t,1),1:nx,t)'...
                *autoCov1Residuals(select,selectLag1)*dg_dx(1:nyIndex(t-1,1),1:nx,t-1)*invALag1;
            % cov(u_t,u_t+1)
            autoCov01Factors(:,:,t) = 1/nyIndex(t,1)*1/nyIndex(t-1,1)*invALag1*dg_dx(1:nyIndex(t-1,1),1:nx,t-1)'...
                *autoCov1Residuals(selectLag1,select)*dg_dx(1:nyIndex(t,1),1:nx,t)*invA;
        end
        
        % Update of nyCount and invA
        invALag1 = invA;
    end
end

% Saving output
out.VarRobust            = VarRobustTheta1;
out.VarHomo              = VarHomoTheta1;
out.VarHomoPooled        = VarHomoPooledTheta1;
out.VarFactorsRobust     = VarFactorsRobust;
out.VarFactorsHomo       = VarFactorsHomo;
out.VarFactorsHomoPooled = VarFactorsHomoPooled;
out.autoCov10Factors     = autoCov10Factors;
out.autoCov01Factors     = autoCov01Factors;
out.autoCov1Residuals    = autoCov1Residuals;
out.varResiduals         = varResiduals;
out.invA_Theta1          = invA_Theta1;
out.A_Theta1             = A_Theta1;
out.PSI                  = PSI;
out.RvHat                = RvHat;
out.RvHat_se             = RvHat_se;
out.RvHatPooled          = RvHatPooled;
out.RvHatPooled_se       = RvHatPooled_se;
out.names                = setup.selectTheta1;
out.residuals            = output.residuals;
out.numBootStep1And3     = setup.numBootStep1And3;
end

